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OBJECTIVES 

The  primary  objectives  of  this  research  is  to  develop  reliable  parallel  algorithms  and  soft¬ 
ware  tools  for  building  improved  reduced-order  models  with  an  emphasis  on  fluid  systems. 

ACCOMPLISHMENTS 


Overview 

Our  research  program  has  focused  on  reduced-order  models  for  large-scale  systems 
with  an  emphasis  on  fluids  and  control  problems.  A  major  accomplishment  has  been  to 
develop  and  analyze  an  algorithm  suitable  for  computing  proper  orthogonal  decomposition 
(POD)  basis  functions  from  large-scale  CFD  data.  The  main  feature  of  this  algorithm  is 
that  it  can  be  used  with  highly  scalable  CFD  algorithms  that  utilize  domain  decomposition 
to  distribute  data  across  compute  nodes.  In  other  words,  it  works  with  data  sets,  allowing 
POD  to  be  used  with  complex  3D  flows.  This  algorithm  was  tested  by  computing  reduced- 
order  models  for  unsteady  flow  past  a  3D  cylinder  computed  on  System  X  (Virginia  Tech’s 
high  performance  supercomputer  which  was  ranked  number  3  in  the  world  in  November 
2003  and  was  still  in  the  top  50  in  2006).  This  computation  would  not  have  been  possible 
if  the  data  had  to  be  collected  on  one  compute  node. 

We  have  also  investigated  a  number  of  traditional  methods  for  reducing  the  computa¬ 
tional  model  in  complex  flows.  These  include  a  multilevel  approach  that  can  be  used  to 
split  the  computation  between  coarse  and  fine  grids.  An  added  feature  of  this  multilevel 
approach  is  that  it  gave  the  first  theoretical  insight  into  the  relationship  between  the  fil¬ 
ter  radius  in  an  LES  model  and  the  required  discretization  size.  Traditionally,  using  the 
relationship  of  two  discretization  points  in  each  direction  inside  a  filter  radius  had  been 
used  based  on  computational  experience.  We  have  also  studied  accuracy  issues  in  comput¬ 
ing  sensitivity  analysis  for  fluids.  This  research  will  be  applied  in  the  follow  on  proposed 
research  incorporating  sensitivity  analysis  with  POD  in  order  to  extend  the  accuracy  of 
reduced  order  models  in  parameter  space. 

For  large  linear  systems,  alternative  projection-based  methods  based  on  a  rational  Krylov 
projection  framework  were  proposed  and  studied.  For  control  problems,  or  model  reduction 
of  linear  input-output  systems,  these  projection  methods  have  optimal  TL2  properties.  This 
projection  framework  has  also  been  applied  to  the  Fourier  model  reduction  methodology 
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as  well.  Toward  the  end  of  this  research  study,  we  investigated  an  interpolatory  framework 
that  provides  insight  on  the  optimal  parameter  values  in  which  to  sample  data. 

We  will  provide  detail  on  the  high-performance  computational  algorithm  for  parallel 
computation  of  POD  basis  vectors  as  well  as  results  from  its  application  to  3D  flow  past  a 
cylinder  at  a  Reynolds  number  of  525. 


The  Filtered  Subspace  Iteration  Algorithm 

An  enabling  method  for  simulating  complex  PDE  dynamical  systems  is  domain  de¬ 
composition  (see  Widlund  and  Tosseli  for  a  good  overview).  The  general  philosophy  is 
to  partition  the  spatial  domain  into  subdomains  and  integrate  the  PDE  in  parallel  on  each 
subdomain.  “Information”  at  the  boundaries  of  these  subdomains  are  communicated  to 
their  neighbors  to  facilitate  the  integration.  There  are  many  approaches  for  implementing 
domain  decomposition  that  we  will  not  discuss  here.  In  practice,  these  subdomains  are 
constructed  so  that  the  size  of  the  discrete  systems  are  nearly  the  same  (for  good  load  bal¬ 
ancing)  and  minimize  the  discretization  size  along  the  interfaces  (for  low  communication 
requirements). 

At  the  end  of  a  simulation,  the  snapshot  matrix  is  distributed  across  several  processors. 
We  present  an  algorithm  for  computing  the  SVD  that  takes  advantage  of  this  data  structure. 
This  is  of  interest  since  it  may  not  be  possible  to  assemble  the  entire  snapshot  matrix  on  a 
single  processor. 

We  assume  that  the  POD  snapshot  matrix  can  be  represented  as 
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The  algorithm  provides  an  efficient  method  for  computing  the  dominant  q  POD  basis 
vectors  where  we  assume 

Tip 

q  «  p  «  Ni  «  N  = 

i= 1 

we  use  caligraphic  notation  {y,U,  V)  to  denote  global  quantities,  subscripts  to  denote  pro¬ 
cessor  number  T),  Y2,  . . . ,  superscripts  to  denote  iteration  number  and'  to  denote  an  ag¬ 
gregated  matrix.  The  algorithm,  known  as  filtered  subspace  iteration  consists  of  two  parts. 
The  first  constructs  a  good  initial  guess  and  the  second  performs  iterative  corrections  to  it. 


Algorithm  1  Initialization 
for  i  —  1:7 
processors  (or  to  processor  0)  end 
Accumulate  V  = 


!0)  =  q  dominant  right  singular  vectors  of  Y)  send  to  other 


vpvp • . .  v<°> 


V!0)  =  q  dominant  left  singular  vectors  ofV. 
broadcast  Vl  0)  to  all  processors  (if  not  locally  computed ) 

Compute  Cil])  =  yVA])  (in  parallel) 

Compute  UA])  =  q  dominant  left  singular  vectors  ofC({>]  using  Algorithm  2  or  3  below. 


This  algorithm  relies  on  the  fact  that  the  snapshot  matrix  comes  from  a  dynamical 
system  and  that  solutions  would  have  the  “separation  of  variables”  structure 

OO 

y(t)  =  te(ta,tb). 

i— 1 

In  other  words,  we  assume  that  the  dominant  time  trajectory  a i(-)  is  prominent  for  each 
portion  of  y\.  Numerical  experiments  indicate  that  there  is  a  strong  correlation  in  the  right 
singular  vectors  of  Yj.  Thus,  this  initialization  provides  a  good  starting  point  for  the  filtered 
subspace  iteration  algorithm  (Algorithm  4). 

One  implementation  note  is  that  the  accumulation  of  V  could  occur  on  every  processor 
(through  broadcasts  of  the  Vy0) ’s).  Hence,  each  processor  would  find  V^01  and  avoid  the 
subsequent  broadcast.  This  avoids  two  communication  steps  at  the  expense  of  one  more 
network  intensive  step.  This  modification  has  not  yet  been  tested.  A  second  note  is  that 
the  Uiiy]  only  needs  to  be  an  orthonormal  basis  for  C,u>  in  practice,  thus  could  be  found  by 
Gram-Schmidt  or  something  equivalent.  However,  using  an  SVD  at  this  steps  gives  good 
agreement  with  the  true  POD  basis  and  hence  is  useful  in  comparisons. 

We  now  present  two  approaches  for  computing  U^\  the  q  dominant  singular  vectors  of 
(an  N  x  q  matrix). 

Algorithm  2  Real  Schur 

1.  Compute  (C^  C^\  requires  one  time  communication  of  approximately  N  x  q. 

2.  Perform  a  real  Schur  factorization  of  this  q  x  q  matrix  product 

(c(j))T  C(i)  =  ZS2ZT. 

3.  Broadcast  Z  and  the  diagonal  elements  of  S  to  all  processors. 

4.  The  multiplication  fZS -1  places  the  correct  rows  qf  U<j  !  onto  the  appropriate  pro¬ 
cessor. 


Note  that  the  product  CTC  above  can  be  performed  in  parallel  with  relatively  little  commu¬ 
nication.  If  C  has  the  following  structure 


CT  = 


then  CTC  can  be  expressed  as  the  sum 

C[c\  +  cZc2  +  ---  +  cZicnp 

where  the  ?'th  term  above  can  be  computed  locally  on  the  ?th  processor.  Thus,  there  only 
needs  to  be  parallel  communication  of  the  np  distinct  q  x  q  matrices. 


Algorithm  3  Gram-Schmidt  Orthogonalize  the  columns  ofC{:,K 


The  result  of  the  initialization  algorithm  is  a  good  approximation  to  the  dominant  q 
dimensional  subspace  for  the  range  of  y.  In  practice,  this  provides  a  great  start  for  the 
iterative  corrections  given  in  the  filtered  subspace  iteration  algorithm  below. 

Algorithm  4  (Filtered  Subspace  Iteration)  for  j=l:Jmax 


1.  Calculate  =  (iT1'3  ^  y.  This  q  x  p  matrix  can  be  accumulated 
processor  or  broadcast  to  every  processor. 


onto  one 


2.  V(j>  =  q  dominant  right  singular  vectors  of  'Dtj).  Possibly  communicate  to  every 
processor. 

3.  Calculate  =  yV<3>  (an  N  x  q  matrix). 

4.  Compute  Uij  ]  as  the  q  dominant  left  singular  vectors  of  C(j>  using  one  of  the  two 
algorithms  above. 

Test  for  convergence. 

Analysis  of  the  filtered  subspace  iteration 

At  the  end  of  the  first  step  in  Algorithm  4,  we  have 

v(3)  =  (uw)Ty 

(since  are  the  dominant  right  singular  vectors  of  D<J>). 

Rearranging  terms,  we  have 


( Uu 


-i) 
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V(j)v(i)l  =  U{j)  E(j) 


By  filling  out  the  orthogonal  complement  to  IT  ^  v>,  we  have 
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Note  that  the  (1,2)  block  above  (in  red)  is  a  zero  block  since  UTy  is  a  rank  q  matrix  and 
the  columns  q  +  1  — >  p  of  are  zero. 

Thus,  we  observe  that 
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Likewise,  after  the  second  step  of  the  filtered  subspace  algorithm,  we  have 
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Again,  the  (1,2)  block  above  (in  red)  is  a  zero  block  since  is  a  rank  q  matrix 

and  columns  q  +  1  — >  p  are  zero. 

Application  of  the  filtered  subspace  iteration 

This  algorithm  was  tested  on  the  three  dimensional  flow  past  a  circular  cylinder  with  a 
fine  mesh  in  the  streamwise  and  crossflow  plane  (128  x  192  x  32)  leading  to  3.14  x  106 
degrees  of  freedom  partitioned  over  128  processors.  Data  was  generated  over  80  time 
snapshots  and  the  above  algorithm  was  applied.  The  first  two  POD  modes  for  the  velocity 
are  depicted  in  Figures  1,  2  and  3.  We  note  that  a  study  of  the  accuracy  of  the  algorithm 
was  performed  in  the  paper  A  Domain  Decomposition  Approach  to  POD  by  the  co-PIs. 


Figure  1 :  Streamwise  POD  Modes  1  and  2 


Figure  2:  Crossflow  POD  Modes  1  and  2 


Figure  3:  Spanwise  POD  Modes  1  and  2 


Continuing  Research 


Two  of  the  co-PIs,  Jeff  Borggaard  and  Traian  Iliescu  along  with  two  postdoctoral  asso¬ 
ciates  (Andrew  Duggleby  in  Mechanical  Engineering  and  Alexander  Hay  in  Mathematics) 
are  using  turbulent  flow  data  to  generate  a  reduced-order  basis  for  improved  reduced-order 
models  for  turbulent  flows.  We  are  currently  validating  the  reduced-order  models  we  pro¬ 
posed  in  challenging  three-dimensional  simulations  of  turbulent  pipe  flows.  We  have  also 
proposed  alternative  approaches  to  the  usual  reduced-order  modeling  approaches,  by  us¬ 
ing  the  Variational  Multiscale  method  and  the  Dynamic  procedure  (based  on  a  two-scale 
computation  of  the  turbulent  model  parameters). 
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Presentations  at  meetings,  conferences  or  seminars 

During  this  project,  we  have  given  more  than  fifty  two  presentations  at  meetings,  confer¬ 
ences  or  seminars.  This  includes  two  tutorial  talks  (a  short  course  and  a  plenary  overview) 
on  reduced-order  modeling.  In  addition,  graduate  students  attended  regional  meetings  and 
gave  presentations  on  their  research  involving  reduced-order  modeling  for  control  prob¬ 
lems. 

C.  Beattie 

1.  SIAM  Annual  Meeting,  New  Orleans,  LA  (July  2005). 

2.  Sandia  National  Lab,  Albuquerque,  NM  (July  2005). 

3.  Ninth  Copper  Mountain  Conference  on  Iterative  Methods  (April  2006). 

4.  Plenary  talk:  Sixth  International  Workshop  on  Accurate  Solution  of  Eigenvalue 
Problems  (IWASEP6),  State  College,  PA  (May  2006). 

5.  SIAM  Annual  Meeting,  Boston,  MA  (July  2006). 

6.  GAMM-SIAM  Conference  on  Applied  Linear  Algebra,  Dusseldorf,  Germany  (July 
2006). 

7.  Workshop  on  Structured  Perturbations,  and  Distance  Problems  in  Matrix  Computa¬ 
tions  (March  2007). 

8.  6th  International  Congress  on  Industrial  and  Applied  Mathematics,  Zurich,  Switzer¬ 
land  (July  2007). 

J.  Borggaard 

1.  Short  Course:  Reduced-Order  Model  Development  and  Control  Design  (with  K. 
Kunisch),  SIAM  Meeting  on  Control  and  its  Applications,  New  Orleans,  LA  (July 
2005). 

2.  SIAM  Annual  Meeting,  New  Orleans,  LA  (July  2005). 

3.  Plenary  talk:  Modeling  and  Optimization:  Theory  and  Applications,  Windsor,  On¬ 
tario  (July  2005). 

4.  ALOSR  Computational  Mathematics  Program  Review,  Long  Beach,  CA  (August 
2005). 

5.  Workshop  on  Large-Scale  Robust  Optimization,  Santa  Le,  NM  (August  2005). 

6.  Austrian  Mathematical  Society,  Klagenfurt,  Austria  (September  2005). 


7.  Florida  State  University,  Tallahassee,  FL  (November  2005). 

8.  Florida  State  University,  School  of  Computational  Science,  Tallahassee,  FL  (January 
2006). 

9.  6th  International  Conference  on  Cooperative  Control  and  Optimization,  Gainsville, 
FL  (February  2006). 

10.  30th  Annual  Conference  of  the  South  African  Society  for  Numerical  and  Applied 
Mathematics,  Stellenbosch,  South  Africa  (April  2006). 

11.  SIAM  Annual  Meeting,  Boston,  MA  (July  2006). 

12.  27th  Annual  Southeastern- Atlantic  Regional  Conference  on  Differential  Equations, 
Greensboro,  NC  (October  2006). 

13.  7th  World  Congress  on  Computational  Mechanics,  Los  Angeles,  CA  (July  2006). 

14.  AFOSR  Joint  Program  Review,  Atlanta,  GA  (August  2006). 

15.  SIAM  Conference  on  Computational  Science  and  Engineering,  Costa  Mesa,  CA 
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16.  14th  International  Conference  on  Finite  Elements  in  Flow  Problems,  Santa  Fe,  NM 
(March  2007). 

17.  Computational  and  Mathematical  Methods  in  Science  and  Engineering,  Chicago,  IF 
(June  2007). 

18.  SIAM  Control  Conference,  San  Francisco,  CA  (June  2007). 

19.  Air  Force  Research  Faboratory,  Wright-Patterson  Air  Force  Base,  OH  (July  2007). 

20.  IFIP  TC7  Conference  on  System  Modelling  and  Optimization,  Krakow,  Poland  (July 
2007). 

21.  AFOSR  Joint  Program  Review,  Fong  Beach,  CA  (August  2007). 

22.  Model  Order  Reduction  Seminar,  MIT,  Cambridge,  MA  (October  2007). 

23.  Plenary  Talk:  28th  Annual  Southeastern- Atlantic  Regional  Conference  on  Differen¬ 
tial  Equations,  Murray,  KY  (October  2007). 

24.  Short  Course:  Farge  Scale  Optimization  and  Design,  DoD  High  Performance  Com¬ 
puting  Program  Office,  University  of  Tennessee  Space  Institute,  Arnold  AFB,  TN 
(February  2008). 

S.  Gugercin 


1.  SIAM  Annual  Meeting,  New  Orleans,  FA  (July  2005). 


2.  Bilkent  University,  Ankara,  Turkey  (August  2005). 

3.  Workshop  on  Large-scale  Robust  Optimization,  Santa  Fe,  NM  (August  2005). 

4.  Massachusetts  Institute  of  Technology  (MIT),  Computational  Prototyping  Group, 
Cambridge,  MA  (October  2005). 

5.  Massachusetts  Institute  of  Technology  (MIT),  Aerospace  Computational  Design  Lab, 
Cambridge,  MA  (October  2005). 

6.  Ninth  Copper  Mountain  Conference  on  Iterative  Methods,  Copper  Mountain,  CO 
(April  2006). 

7.  Conference  on  Adaptive  Model  Reduction  Methods  for  PDE  Constrained  Optimiza¬ 
tion,  Rice  University,  Houston,  TX  (May  2006). 

8.  Joint  GAMM-SIAM  Conference  on  Applied  Linear  Algebra,  Duesseldorf,  Germany 
(July  2006). 

9.  SIAM  Conference  on  Computational  Science  and  Engineering,  Costa  Mesa,  CA 
(February  2007). 

10.  Workshop  on  Structured  Perturbations,  and  Distance  Problems  in  Matrix  Computa¬ 
tions,  Bedlewo,  Poland  (March  2007). 

11.  Koc  Univeristy,  Turkey  (July  2007). 

12.  Sixth  International  Congress  on  Industrial  and  Applied  Mathematics,  Zurich,  Switzer¬ 
land  (July  2007). 

T.  Iliescu 

1.  The  58th  Annual  Meeting  of  the  Division  of  Fluid  Dynamics  of  the  American  Phys¬ 
ical  Society,  Chicago  (November  2005). 

2.  Workshop  on  Collaborations  in  the  Mathematical  Geosciences,  poster  presentation, 
Research  Triangle  Park,  NC  (October  2005). 

3.  Southeastern  Section  MAA  and  SIAM  Southeast  Atlantic  Section  Joint  Meeting, 
Auburn,  AL  (March  2006). 

4.  University  of  Ottawa,  Department  of  Mathematics  and  Statistics,  Ottawa,  Canada 
(May  2006). 

5.  Universita  di  Pisa,  Department  of  Applied  Mathematics,  Pisa,  Italy  (May  2006). 

6.  Laboratory  for  Modeling  and  Scientific  Computing  MOX,  Department  of  Mathemat¬ 
ics  ”F.  Brioschi”,  Politecnico  di  Milano,  Milano,  Italy  (May  2006). 

7.  SIAM  Conference  on  Computational  Science  and  Engineering,  Costa  Mesa,  CA 
(February  2007). 


8.  6th  International  Congress  on  Industrial  and  Applied  Mathematics,  Zurich,  Switzer¬ 
land  (July  2007). 

9.  Plenary  Talk:  Sandia  CSRI  Workshop  on  Mathematical  Methods  for  Verification 
and  Validation  (August  2007). 

10.  AMS  Spring  Central  Meeting,  Bloomington,  IN  (April  2008). 

Interactions  and  Transitions 

Air  Force  Research  Laboratory ,  Wright-Patterson  Air  Force  Base,  OH 

Jeff  Borggaard  spent  the  summer  at  AFRL/VACA  working  with  Chris  Camphouse  and 
James  Myatt  on  a  flow  control  problem.  Efficient  POD  software  using  algorithms  devel¬ 
oped  in  this  research  and  applicable  to  practical  engineering  reduced-order  model-based 
control  algorithms  was  provided  to  the  AFRL  researchers. 

Daniel  Sutton,  Masters  student  in  Mechanical  Engineering,  and  partially  supported  on  this 
project,  spent  a  summer  at  the  Air  Vehicles  Directorate  working  with  Siva  Banda,  Chris 
Camphouse  and  James  Myatt:  He  worked  on  computing  low-order  models  based  on  Proper 
Orthogonal  Decomposition.  He  studied  the  use  of  POD  for  coupled  systems.  The  PI  and 
Lizette  Zietsman  made  a  two  day  visit  to  the  lab  during  Daniel’s  internship  to  discuss  his 
research. 

[Contacts:  Chris  Camphouse  (937)  255-6326,  James  Myatt  (937)  255-8498] 

Synergistic  Activities 

1.  Serkan  Gugerin,  with  Karen  Willcox  from  MIT,  organized  a  two-part  minisympo¬ 
sium  on  Model  Reduction  at  the  SIAM  Annual  Meeting,  in  Boston,  MA,  July  10-14, 
2006. 

2.  Jeff  Borggaard  and  Traian  Iliescu  co-organized  Emerging  Finite  Element  Methods 
for  Complex  Flow  at  the  2007  SIAM  Conference  on  Computational  Science  and 
Engineering,  Costa  Mesa,  CA  (February  2007). 

3.  Jeff  Borggaard  co-organized  Model  Reduction  Methods  for  Flow  in  Porous  Media 
at  the  2007  SIAM  Conference  on  Mathematical  and  Computational  Issues  in  the 
Geosciences,  Santa  Fe,  NM  (March  2007). 

4.  Joseph  Ball  and  Christopher  A.  Beattie  and  Serkan  Gugercin  from  Virginia  Tech., 
Athanasios  C.  Antoulas  from  Rice  University,  and  Tryphon  T.  Georgiou  from  Univer¬ 
sity  of  Minnesota  are  the  organizers  of  the  18th  International  Symposium  on  Math¬ 
ematical  Theory  of  Networks  and  Systems  (MTNS  2008)  in  Blacksburg,  VA,  July 
2008.  There  are  several  sessions  on  reduced-order  modeling.  Presentations  by  the 
co-PIs  and  graduate  students  on  our  reduced-order  modeling  research  are  scheduled. 

Honors/Awards 

Jeff  Borggaard  was  awarded  an  ASEE  Summer  Faculty  Fellowship  and  spent  May-July 
2007  at  the  AFRL  Control  Sciences  Center  of  Excellence. 

Serkan  Gugercin  was  awarded  an  NSF  Early  Career  Award  in  Computational  and  Applied 
Mathematics,  2007. 


